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Abstract 

The Hohenberg-Kohn theorem and Kohn-Sham procedure are extended to functionals of the 
localized intrinsic density of a self-bound system such as a nucleus. After defining the intrinsic- 
density functional, we modify the usual Kohn-Sham procedure slightly to evaluate the mean-field 
approximation to the functional, and carefully describe the construction of the leading corrections 
for a system of fermions in one dimension with a spin-degeneracy equal to the number of particles 
N . Despite the fact that the corrections are complicated and nonlocal, we are able to construct 
a local Skyrme-like intrinsic-density functional that, while different from the exact functional, 
shares with it a minimum value equal to the exact ground-state energy at the exact ground-state 
intrinsic density, to next-to-leading order in 1/N. We briefly discuss implications for real Skyrme 
functionals. 
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I. INTRODUCTION 



The Hohenberg-Kohn theorem [jj and the associated Kohn-Sham procedure |2j have 
helped to shift the interpretation of Skyrme-mean-field theory in heavy nuclei. Nowadays, 
instead of focusing on a Skyrme nucleon-nucleon interaction, theorists start with a Skyrme 
(or relativistic) mean-field-like density functional, modifying its parameters to fit nuclear 
properties in ways that are not obviously consistent with the existence of an underlying 
effective two-body interaction. The Hohenberg-Kohn theorem is often used to justify this 
step, the idea being that Skyrme functionals are an approximation to a universal Kohn-Sham 
density functional — written in terms of single-particle orbits and leading to Hartree-like 
equations — that is guaranteed to include all correlations. 

Because the nucleus is a self-bound system, however, the usual Hohenberg-Kohn theorem 
is not really a justification. It states that there exists a universal functional of a system's 
laboratory density that has a minimum at the ground state density, and that the value of the 
functional at the minimum is the ground-state energy. The adjective "universal" means that 
the functional can be written as a sum of a part that is independent of whatever external 
one-body potential the system is in — the universal part — and a simple term J V(r)p(r)dr, 
where V is the one-body potential and p is the usual laboratory density. The density of 
interest in nuclear physics, however, is the intrinsic density, i.e. the density relative to the 
nuclear center of mass, not relative to a fixed point in the lab. It is this intrinsic density 
that is localized and measured in scattering experiments. By contrast, the laboratory nuclear 
density of an isolated nucleus is spread out evenly over all space because the nuclear center 
of mass is free to move anywhere. Thus, the universal Hohenberg-Kohn functional has a 
minimum at a constant spread-out density that contains no information about the localized 
intrinsic density. The guaranteed existence of this functional is irrelevant to nuclear physics. 
The same is true of the usual Kohn-Sham procedure, which represents the density as a 
sum of squares of single-particle wave functions, leading to a Hartree-like equations for the 
ground-state density and energy. The single-particle orbitals for the ground-state laboratory 
density are plane waves. 

Because the ground state wave function of the nucleus factors into center-of-mass and 
intrinsic parts, it is possible - - by removing the center-of-mass kinetic energy from the 
Hamiltonian and adding a potential such as Vqm = (2j r i) 2 that * s minimized when the 
center of mass is at the origin — to construct a related Hamiltonian for which the ground- 
state laboratory density is the same as the intrinsic density. One could then try to produce 
an ordinary density functional for this new Hamiltonian, assured by the Hohenberg-Kohn 
theorem that its minimum would yield the true intrinsic density. By making the potential 
very weak, one could even make the functional yield the true ground-state energy when 
minimized. This procedure, however, would single out one point in space and would therefore 
look very different from Skyrme mean-field theory, in which the functional is minimized by a 
family of Slater determinants, related to one another by translation. Furthermore, the local- 
density approximation (LDA) and extensions, through which Skyrme functionals might be 
derived j.3|, could easily destroy the conditions that force the intrinsic density to be the same 
as the laboratory density. The whole Skyrme enterprise would make much more sense if we 
knew that we could construct an exact Kohn-Sham functional of of the intrinsic density 
rather than the laboratory density, without modifying the Hamiltonian. 

These considerations, and related ones connected with deformation and rotation, lead us 
to ask whether analogs of the Hohenberg-Kohn theorem and Kohn-Sham procedure exist 
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for the intrinsic density. Rotation is more complicated than translation because it doesn't 
decouple from internal motion and because the intrinsic rotational density is not directly 
observable, so for now we worry only about translation and define the intrinsic density 
simply as the density with respect to the center of mass. The operator corresponding to 
that quantity is 

Pi{r) = p(r + r CM ) , (1) 

where p is the usual laboratory density operator, the center-of-mass position operator f cm 
is defined by 

1 - 

*™=jvX) f< ' (2) 

t=i 

and the hats over pj, p, r i} and Tqm are to emphasize that the quantities are operators, 
unlike the variable r in Eq. The intrinsic density, unlike its lab-frame counterpart, is 
an iV-body operator, complicating density-functional theory for self-bound systems. 

In the next section we address the question of whether there is an analog of the Hohenberg- 
Kohn theorem guaranteeing the existence of an intrinsic "universal" intrinsic-density func- 
tional with a minimum at the exact ground-state intrinsic density and energy. Section 11111 
argues that a Kohn-Sham-like procedure can be used to find the exact intrinsic density, and 
outlines the construction of the corresponding Kohn-Sham equations. We will work with 
a specific one-dimensional model to simplify the formalism, but will carry the construction 
as far as possible; as a result this section is the longest and most technical. In section IIVI 
we produce a phenomenological Skyrme-like functional that, while not exact away from the 
minimum, does yield the model's exact intrinsic density and energy at its minimum to good 
accuracy. We then briefly discuss the implications of these results for real three-dimensional 
nuclear physics. 



II. HOHENBERG-KOHN THEOREM FOR INTRINSIC DENSITY 



There is indeed a Hohenberg-Kohn theorem for the intrinsic density. In fact, Valiev and 
Fernando [H have shown — assuming a nonrelativistic Hamiltonian, and using an effective- 
action formalism invented in Ref. - that given any Hermitian operator Q(r), one can 
construct an energy functional E[Q] that is universal in a sense to be discussed and has a 
unique (and correct) minimum at Q(r) = (Q(r)) gs , the ground-state expectation value of 
Q(r). Introducing a "source" J(r) that couples to Q(r) and alters the ground-state energy: 

£[J] = (H + J J(r)Q(r)dr) J gs , (3) 

one can define the functional E[Q] as the functional Legendre transform with respect to J. 
In other words 

E[Q]=£[J[Q]]-J[Q]oQ , (4) 

where J[Q]{r) is the particular source function that makes {Q{r)) J ?s = Q(r), and we have 
used the convention A o B = f A(r)B(r)dr. It's not hard to show [6( (and this is an alter- 
native and perhaps clearer definition) that E[Q] is the minimum of (H) over all normalized 
iV-particle wave functions that have (Q(r)) = Q(r). Since (H) is smallest in the ground 
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state, the functional is minimized when Q = (Q) gs - It is universal in that the addition of an 
arbitrary term V o Q to the Hamiltonian merely adds V o Q to E[Q]. 

All this means, in particular, that for nuclei there is a universal "intrinsic-density func- 
tional" E[pi] with a minimum at the ground-state intrinsic density, and that we know in 
principal how to construct it. It's worth noting that the functional is not universal under 
the addition of an ordinary one-body potential ^ V(^i) to the Hamiltonian, but rather 
the addition of a generically iV-body potential £\ V{fi — tcm) that affects only intrinsic 
structure. 

The important question, though, is how to incorporate orbitals into this framework in a 
way that leads to a constructive procedure for the energy functional. In other words, can 
one derive a Kohn- Sham-like (e.g. Skyrme-like) functional for the intrinsic density, even in 
principle? In the next section we show that one can. 



III. ORBITALS AND KOHN-SHAM-LIKE EQUATIONS 
A. Framework 



Refs. |4j, [7[ introduce Kohn-Sham orbitals in functionals of the ordinary density (in sys- 
tems that are not self bound) through the perturbative "inversion method" In this 
approach one divides J[p], where the tilde here and below distinguishes the quantity from 
the related one to be used here for the intrinsic density pi, into a piece Jo[p], the Kohn- 
Sham potential that forces a noninteracting system to have the density p, and the rest. The 
entire functional can be constructed from J an d expressed in terms of the noninteracting 
orbitals <pk and energies that Jo determines. In our notation, under the assumption that 
S[J) can be divided into successively smaller contributions £[J] = S [J] + E\ [ J] + ■ ■ ■ , the 
energy-density functional E[p\, has the form 0,0 

E[p}=E [p}+E mt [p} , (5) 

where 

E [p] = S [J ] -J op, E int [p] =S 1 [J } + -- - , (6) 

and we have truncated -E\ nt [/>] at lowest order. To find the ground-state density and energy 
one rewrites the minimization condition SE[p]/5p = 0, as an implicit Kohn-Sham equation 
of the form 

^ = JM , (7) 

bp 

where we have used Eq. (0) and the fact that 5E Q /5[p] = — J , which follows from the defini- 
tion of E Q as a Legendre transform. The Kohn-Sham procedure then consists of choosing a 
trial potential J , solving the single-particle Schroedinger equation to obtain the correspond- 
ing (pk and e/c, computing E int and 5E int /Sp in terms of these quantities, and then iterating 
until the energy and density converge. The procedure amounts to solving Hartree equations 
with the mean- field potential V[p] = 6E int /Sp. 

We can apply this procedure to our intrinsic functional E[pj], even though pj is an iV- 
body operator, provided we start at a slightly different point: mean- field theory. To make 
the discussion less abstract, we refer explicitly to a simple one-dimensional system. Before 
constructing an energy functional, we describe the system, its exact ground-state energy 
and intrinsic density, and the mean-field approximations to these quantities. 
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B. Model in one dimension 



We consider N fermions with spin-degeneracy N interacting via an attractive two-body 
delta function (h— 1): 

1 N d 2 

i=l 1 i<j 

This model has the nice feature of being exactly solvable. The ground state is bound and 
the corresponding energy and intrinsic density are known jol llfll]: 

= ^mg 2 (N 3 - N) (9) 




n(N\) 2 e~ nNm9 ^ 



(N + n-l)\(N -n-l)\ 
3 



2cosh 2 (^p; 




mg\x\ 



> 



(10) 

3 IniV 
N 2 ' 



The ground state differs from that of a nucleus, however, in that the system shrinks to 
a point as N increases and the binding energy goes like N 3 ; there is nothing resembling 
saturation. A local-density approximation for such a system, which has no uniform limit, is 
hard to imagine. 

Ref. [Icj analyzed the Hartree-Fock approximation, to which we refer as the mean-field 
approximation, for the ground state. The approximation is improved when the center-of- 
mass kinetic energy is subtracted from the Hamiltonian, as e.g. in the no-core shell model 
[ill ]. Subtraction makes a localized ground-state possible, since the center-of-mass is no 
longer required to be in a plane-wave state. Because the spin degeneracy is equal to the 
number of particles, whether or not one subtracts, the spatial mean-field orbital <p(x) is the 
same for every particle. The mean-field equation for 0, with no external source and the 
center-of-mass kinetic energy subtracted, is 

A " V(z) ~{N- l)g\<P(x)\ 2 ^x) = e<f>{x) (11) 



2mN 



Because of translational invariance, the solutions make up a degenerate set with each member 
centered at a different point x = a in space: 



The average position x, given by 

*=/x|#M|>«fa, (13) 
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is just x = a for the mean- field solutions. The mean- field energy is then 



mf 



N(N - 1) J dx 
N(N-l)g 



^(*)l a -fl*(*)l 4 



Ne + 



24 



(x)\ 4 dx 



-N 2 (N-1) 



(14) 



Interestingly, the mean-field energy is correct at leading order in the small quantity 1/N, 
that is, the 0(N 3 ) term is right. At 0(N 2 ) there is an error; the exact energy in Eq. (Julias 
no term of that order. The mean-field Q(N 2 ) energy is reduced by the subtraction of the 
center-of-mass kinetic energy, but, as Ref. [l2[ shows, one needs to include ring diagrams 1 
(i.e. the RPA correlation energy) to cancel it completely. 

The mean-field laboratory density N(j) 2 (x) is also correct in some sense in leading order; 
because center-of-mass motion becomes irrelevant at large N, the mean-field laboratory 
density reproduces the exact intrinsic density there, assuming the center of mass is at the 
origin The mean-field intrinsic density can be computed by expanding the operator 

x) around x + x in powers of xcm — %, leading, after some algebra, 



= p(x + X + X C M 

to 



(Pi(*)) 



mf 



c N N 



(cnx + x)\ 2 + J^^l^x + x)\ 2 + <D(± 



N 



\<t>{x + x)\ 2 + - 



(x + x)\ 2 + £-7~l 

ax 



(x + x)\ 2 



(15) 



with 



x — X) 



'■\(f)(x)\ 2 dx 



N 



cn 



N-l 



(16) 



Beyond leading order, this expression gives a better approximation to the exact intrinsic 



density than does the mean- field laboratory density N\ 



as can be seen, for example, 



by using it to compute the (x 2 )i ntr . The result, (x 2 )i ntr = /x(l — jt), is the same as the exact 
one. 

With these results in hand, we can move to the construction of a Kohn-Sham-like intrinsic- 
density functional for the model. 



C. Intrinsic Density Functional at Mean-Field Level 

The problem with applying the Kohn-Sham procedure directly to the intrinsic density is 
that the source term, Jo o pj, is an iV-body operator. We cannot treat the source exactly 
and still have a noninteracting system, the starting point for the inversion method discussed 



One might mistakenly infer from Ref. |12( that ring diagrams with more than two rings contribute very 
little to the 0(N 2 ) correlation energy. That can appear to be so, but only if one modifies the Hamiltonian 
(by subtracting the center-of-mass energy) in the Hartree-Fock calculation but neglects to do so in the 
RPA. 
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above. Instead, we start with the mean-field approximation, not just for the true ground- 
state, but also for the ground state with the addition of the term Jq o pj. In other words, 
Jo[pi] is now the source function that forces the intrinsic density to be pi when the system 
is treated in mean-field approximation. The resulting mean-field single-particle orbitals will 
play the same role in obtaining the exact intrinsic functional as Kohn-Sham orbitals for a 
noninteracting system do in the usual approach, and the inversion method will go through 
without any other modifications. Of course, to apply it we must assume that corrections 
to mean-field theory are perturbative. For the model just described that is certainly the 
case, as we shall see. For real Skyrme functionals, it is not unreasonable to assume the 
existence of an effective Hamiltonian for which Hartree-Fock is a good starting point, at 
least in spherical nuclei. 

Constructing the mean- field intrinsic-density functional as the first term E [pi] in the 
intrinsic version of Eq. (jSJ) can be done without solving any equations. We are supposed to 
start with 



£ [J ] = N(N-l) / dx 



2mN' v " 2 



+ ^o° (P/)mf , (17) 



where (/3/) m f is given by Eq. (fT5|) and we have taken the useful but not essential step of 
subtracting the center-of-mass kinetic energy (though without explicitly writing a more 
complicated two-body center-of-mass term that has no effect on bound states) to improve 
the mean-field approximation. But since J is supposed to be fixed so that (p/) m f = Pi and 
then Jq o p l subtracted from £o[^o] ; the mean-field functional looks just like the first line of 
Eq. (HI: 



E [ Pl ]=N(N-l) j dx 



' : |0'(x)| 2 -g|0(x)| 4 



2mN' v yi 2' 



(18) 



This simple Skyrme-like object, when written in terms of the orbital 0, is the same as 
the mean-field functional one would get for the ordinary one-body density, but now is a 
functional of pi (up to the location of the center of mass x, which must also be specified) 
rather than p. The relation between the two functionals is not surprising because at mean- 
field level the Slater determinant constructed from orbital is the system's wave function, 
and the energy functional — the energy corresponding to that wave function — is just the 
expectation value of the Hamiltonian, no matter what observable we write the functional 
(and 0) in terms of. (When we go beyond mean-field theory, center-of-mass motion makes 
the functionals of the lab and intrinsic density quite different.) 

In more detail, is a functional of pi in the following sense: pi determines Jq through 
the solution to 8£q/5Jq = pi (Ref. Q shows that the solution exists and is unique) and 
Jo determines 0, again up to its overall location, through the solution to the mean-field 
equations, which we haven't yet written, that come from varying £o[Jo] with respect to 
(with a constraint on its norm). With and Eq functionals of pi, one can then reproduce the 
mean- field ground-state energy and intrinsic density by solving the equation 5E /5pi = 0. 
Because this equation says that Jq = 0, finding the minimum corresponds to solving the 
source-free mean- field equation, i.e. the Hartree equation Eq. (|11|). and then computing pi 
from Eq. (|15|). This Hartree equation is the leading approximation to the Kohn-Sham-like 
equation we seek. 

All this may seem pedantic, but it contains an important point about the relation of to 
pi: the condition 5£q/5Jq = pi simply states that pi will be given in terms of by Eq. (fT5j) . 
with pi replacing (p/) m f. Thus, though pi(x) is just iV|0(a;-|-x)| 2 to leading order in 1/N, it 
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has higher-order corrections. This fact is apparently important for real Skyrme functionals, 
which fit data best when a prescription |l3( for pi that approximates the real-world analog 



of the Eq. (fT3j) is used to calculate observables such as the rms radius [14(. Using the 
uncorrected expression pi(x) = |0jt(x + x)| 2 for this purpose produces noticeably poorer 
results. 

It actually is possible to formulate an intrinsic-density functional theory in which pi = 
l^kix + x)\ 2 without corrections. In the simple model we can do so by neglecting the 
0{1/N) terms in Eq. (fTo]) when defining £o[Jo], and taking their effects into account at the 
next level of approximation. Such a step would in fact simplify much of the formalism to 
come, e.g., the mean-field equation for the orbital <ft i n the presence of the many-body source 
Jo- And when expressing the leading-order intrinsic-density functional in terms of <fi we'd 
still end up with Eq. (fTHj) . That object, however, would be a different functional of pi than 
before and would produce poorer results for the ground-state intrinsic density. We will not 
elaborate the corrections to the functional necessary to restore the lost accuracy, lest this 
paper get even longer than it is. Instead, using the relation between pi and (f) given by Eq. 
(I15|) . we proceed to take up corrections that are genuinely beyond mean- field theory. Though 
the lack of something like the LDA prevents us from actually calculating the corrections, 
either analytically or numerically, we describe how to do the latter, in some detail. 



D. Beyond Mean-Field Theory 



To go beyond the relatively straightforward mean-field approximation to E [pi] , one would 

need to calculate E- mt [pi] = Ei[pi] + • • • = £i[Jo] + (We content ourselves with E\ 

here, though we could go further.) Then, provided one could calculate SEi/Spi(x), one 
could follow the procedure Kohn-Sham procedure outlined above, rewriting the minimization 
condition 5E[pi]/5pi in the form 



Jnix) 



S£i[Jq_ 
5pi(x) 



(19) 



and solving by iteration, starting with a guess for Jo. We will call Eq. ()19|) the Kohn-Sham 
equation, even though strictly speaking it doesn't involve Kohn-Sham orbitals. 

Given a Jo, one could find Si in the following way: The first step would be to find the 
associated mean-field orbital by solving the mean-field equation, which we finally have to 
write down, for the system in the source J o pj. We obtain the equation by varying in 
Eq. (|T7|) and using Eq. (fTo]) . Since x in the latter is itself a functional of 0, the result is 
complicated: 



<t>"{x) 
2mcN 
/' 



(iv-i)^)! 2 - JS(x) + 



[x — X) 



C N 



J*'(y)my)\ 2 dy 



(p(x) 



(20) 



2N 



[X — (X — X) 



Jr\y)\m\ 2 dy + 



(x — X) 



r "{y)\<i>{y)\ 2 dy 



4>{x) = ecf>(x) 



where Jq(x) = Jo(^f ) and we have been (and will continue to be) cavalier about factors of 
in terms of 0(1/N), since they make a difference only in terms of 0(1/N 2 ), which we have 
dropped. They could be included in a more accurate functional. We have also added x's in 
the fourth and sixth terms to keep translational invariance manifest; they would otherwise 
be absorbed into e. 



S 



Next, one would need an approximation in which to evaluate/define E\. As mentioned 
above, Ref. [lij demonstrates that in the absence of a source the largest corrections to the 
mean-field approximation in a 1/N expansion come from summing ring diagrams. The same 
statement does not apply to the intrinsic density functional at all possible densities, but the 
ring sum, which is the leading correction to the mean-field in a "loop-expansion" of 
S[J], should yield the 1/N corrections near the ground-state intrinsic density and makes a 
natural choice for Si. 

To obtain an expression for the ring sum, we work with the one-body density matrix p a b 
[l6l |. We can use it, for starters, to represent the contents of Eq. ([20)1 . The left-hand side 
of that equation, with a delta function S(x — x') replacing <p(x), gives the coordinate-space 
representation of the mean field h(x, x')[Jo] ■ We can write this h (with spin now included) 
in an arbitrary space+spin basis in terms of p ab : 



h a b[Jo\ 



d£o\M 

dpba 

1" (^0)06 

1 



(21) 



1 



Tr[(J*y P ](x-xU 



2N 2 



c N N 

(»N{Ji): b -^[{4)'"p\ {x 



x) ab + Tr[(J*y P ] [{x 



x" 2 



}ab) 



where Tr represents a trace over matrix elements, and we have only omitted the terms 
that do not depend on Jo in the second and third lines. Hartree-Fock theory corresponds 
to finding the basis that make h diagonal. Again, we have inserted two x's to manifest 
translational invariance; they just add a constant to all eigenvalues, without altering the 
Hartree-Fock basis. 

The ring sum, in these terms, is just the RPA correlation energy calculated with the 
effective two-body interaction [l6[ 



V r , 



ab,cd 



dh ac [J ] d 2 S[J ] 



dpdb dpcadpdb 
In coordinate space, this interaction takes the form 

V(XI,X 2 ) = —S(XI ~X 2 ) + 7T^T7P1P2 

cn 2mN 



(22) 



(23) 



1 

'N 



1 r 



(xi — x) Jq(x2 — x) + (x2 — x) Jq(xi — x) — (xi — x)(x2 — x) ( / Jq(i/ — x)\4> 2 {y)\dy 



where the term with pj = —id/dxi is the two-body center-of-mass Hamiltonian, matrix 
elments of that term alone are to be antisymmetrized, and we have replaced some more 
complicated terms with an ellipsis. The extra x's we added to h (and to Eq. (J23|) ) have no 
effect because the particle-hole RPA involves only terms h ac with a ^ cm. Eq. (|22jl and x is 
just a number, with no off-diagonal matrix elements. 

The sum of all the ring diagrams can be expressed as |l7L 111 



Ei[ Pl \ = Si[J ] = ^- J Re (Tt ln(l - R(iu)V) + R(iuj)V ) du , 



(24) 
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where the coordinate-space matrix elements of R(ui) make up the "unperturbed" response 
function R(xx, x-i, x[, x' 2 ; u>) for the system in the Jo-dependent mean field. (The spin indices, 
which we've omitted, just contribute factors of N in Eq. (|24|h) For zero-range interactions 
and our source we need these matrix elements only at and there 

the response can be written (for real) as: 



R(x, x';uj) = cj){x)(f){x') 



(x\ - \x') + (x'\ \x) 

uj + e + irj — h[J ] ' — io + e + ir] — h[J ] 



(25) 



The matrix elements are just single-particle Green's functions for particles in the mean field. 
Equation ()24|) is then the RPA correction to the intrinsic density functional. 

One advantage of our definition of intrinsic density (for others, see, e.g., Refs. floL 
is that £[J] is translationally invariant; moving 4>(x) doesn't change it because it doesn't 
change pi, the quantity to which J couples. As a result, the RPA equations have a zero 
mode. With the center-of-mass kinetic energy subtracted subtracted from the Hamiltonian, 
the zero mode contributes nothing to Ex, but that convenience is not an essential part of the 
treatment. We could actually start with the full Hamiltonian, with no center-of-mass-energy 
subtraction, and recover the same result for E + E. The mean-field energy E [pi] would 
be different, but the difference would be made up by the contribution of the RPA zero 
mode to Ex, the ring sum in Eq. (j24j). Further corrections might require mixing of mean- 
fields with different values of x to account for center-of-mass motion, but at RPA order a 
single- mean field with a proper treatment of the zero mode is sufficient, even without the 
trick of removing the center-of-mass energy from the Hamiltonian. 

Despite this nice feature, Eq. (J24*Jl is still a complicated implicit representation of Ex[pi]. 
To solve the Kohn-Sham equation, Eq. (|19|) . we need the derivative of this functional. Taking 
the derivative is the most numerically involved step in the entire process 2 . References |4j,|7J 
suggest using the relation 

5Ex[pi] _ 5£x _ 5 J ^ 



o 



5pj (a;) 5J 5pi(x) 

In our model the two terms on the right-hand side can be manageably combined, particularly 
if we're willing to sacrifice the quality of the functional away from the minimum, a step that, 
as we shall see, has no effect on the Kohn-Sham solution. 

To get one can start by writing in terms of pi and x through the perturbative 
inversion of Eq. (JTSj) (assuming <fi real): 



m = (ckN) -»> J P ,( X -^) - £rf(^> + Oi±) , (27) 

Then one could solve Eq. (j2U|) in reverse to obtain J — e in terms of 0. (Jo is needed only up 
to a constant to evaluate higher-order corrections.) Though Eq. (j2~U|) for J is a complicated 
differo-integral equation, we can simplify it without harm. Near the ground-state value of pi, 
the source J , which is close to SEx/Spi there, is smaller by 0(1/ N) than (N — l)g|0| 2 , the 
Jo-independent part of the mean- field equation. Ex is already smaller than Eq by 0(1/ N) 
near the ground state, and the contributions to E\ of terms in which J is multiplied by 



2 Ref. |20j discusses some of the difficulties in "orbital-dependent functionals" for atomic and condensed- 
matter physics. 
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1/N are therefore down from Eq by 0(1/N 2 ). At that order we're already missing terms 
because of the restriction of the relation between pj and <fi in Eq. (j!5|) to 0(1/N). Thus, in 
treating E\ or its derivative we can neglect the entire second line of Eq. ()2U|). (and even the 
last term on the first line because the corresponding to the ground state is symmetric) 
without changing the Kohn-Sham solution at 0(1/N). This converts the differo-integral 
equation for J into a simple algebraic one, yielding 

J °w-^i9r +iV * w|2 - (28) 

With this result and Eq. (}2"Tj) it's easy to get 5J (y)/5pi(x) through the chain rule; it 
would vanish for x ^ y with the approximations just mentioned. One could then extract 
5Ei[pj]/SJ (y) by making local variations in in Jo(y) and calculating the changes in E^pi] = 
£\[Jq] through Eq. (j2I|). One could simplify that equation as well by neglecting all but the 
analogs of the first term after the ellipsis in Eq. (j21)l and the first two lines of the effective 
interaction, Eq. ()22|). without affecting the functional near its minimum at next-to- leading 
order. 

That leaves, finally, the Kohn-Sham equation itself. As mentioned already, one can solve 
it iteratively. For our intrinsic-density functional that means starting with a guess for J , 
solving the mean-field Eq. (}2T)j) to obtain and e, evaluating Ei via Eq. ()24|) and then 
SE 1 /5pi(x) as just discussed, resetting J to equal 5Ei/5pi, and repeating until J , E\, and 
cf) (and therefore pi) converge. This is equivalent to solving the mean-field equation (|2T)|) 
with J replaced by SEi/Spi. And again, because even the leading order terms in Ei are 
down by 0(1/N) from E Q near the minimum, we can throw away the higher-order terms 
in that equation without affecting the relevant part of E\. We can even neglect them in 
constructing E [J [pi]] near the ground state, because in that functional corrections to the 
mean-field energy, Eq. are second order in the small source Jo- One can get the entire 
density functional correct at next-to-leading order near the minimum with a Kohn-Sham 
equation that takes simple Hartree form! 

For real nuclei, in which there is more than one spatial orbit, there are no analogs of 
Eqs. (J27J) and and our procedure for evaluating 6J /5pi won't work. That quantity is 
the inverse of a generalized linear response function for pi at uj = 0, and will not usually 
vanish at y ^ x, as it did in our machinations above. It be evaluated numerically, but at 
a high computational cost. Sub leading terms can still be neglected in this inverse response 
function, however, if the only goal is to include the leading corrections to the mean-field 
functional near the minimum. That would mean ignoring the difference between SJo/Spi(x) 
in Eq. (|26J) and the ordinary RPA inverse res pon se function, which for a given Jo can be 
calculated in a relatively straightforward way |2l|. Though 8£i/8J , the other quantity in 
Eq. (J2HJ), might still not easy, it could be evaluated as already described or (particularly for 
real and more complicated functionals) via the relations in Ref . • 

If one wanted to avoid even that difficulty, one could make still more approximations 
without changing the solution. One could even go so far as to evaluate SEi/SJo and SJo/Spi 
at J = 0, i.e. at pi = (pi(x)) m f = (mgN 2 /4)cosh~ 2 (Nmgx/2) (the mean-field approxima- 
tion to the exact value). That would actually remove Jo from the right hand side of Eq. 
(J2HJ), thus changing the self-consistent procedure represented by Eq. (fHJj) into a one-step 
formula: 



J (x) 



5Jo(y) 



S£i 

x 



Spi(x)) j =o 5J (y) 



dy , (29) 

Jo=0 
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which, after multiplying both sides by the response function Spi(z)/SJo(x) (which equals 
8pi(x)/8Jo(z)), and integrating over x, implies 

P/0) = (pi(x))mf\ Jo=0 + [ T7T~T - , J o(y)dy + 



SJ (y) 



i <- r w ^1 
{Pl{x))mf + 



J =0 



Jo=0 

(30) 



8J (x) 

(The manipulations serve to obtain the second line from the first, which is just a Taylor 
expansion.) Evaluating Eq. (|3T)j) would just be equivalent to calculating the correction to 
the mean-field ground-state intrinsic density by summing rings generated by the attractive 
delta function and the first-order insertion of an effective one-body potential (the third 
term in Eq. (j2TJjO or two-body interaction (the second line of Eq. (J221))) with J(y) replaced 
by 5(x — y) in both cases because of the functional derivative. That is eminently doable 
and would reproduce the second line of Eq. (JTUJ), the expression for the exact ground-state 
intrinsic density. 

In our model, however, no matter where in this chain of approximations we decided 
to stop, the intrinsic-density functional would not be local, i.e., it would not look like a 
Skyrme functional. For that, one would seem to need some form of LDA. The answer to 
the question at the end of the introduction, then, is that we can write the intrinsic density 
functional in terms of orbitals, but it's tough, in this model anyway, to make it Skyrme-like. 
Nonetheless, the successive approximations do bring an important point to light: what we 
really care about is not so much the density functional for arbitrarily stressed nuclei as it is 
the ground-state energy and intrinsic density itself. In fact, since the nucleus is never subject 
to a deforming potential, except in nuclear matter where wave functions are translationally 
invariant, one might argue that that is all we care about. It actually seems a bit strange to 
apply density-functional theory, the main feature of which is "universality", to a problem 
in which there are no external influences to treat universally (though density functionals do 
arise naturally without the accompanying potentials in the density-matrix expansion [3]). 
In any event, a functional that does not exactly match the definition, Eq. (j3J) but does give 
the correct energy and intrinsic density when minimized, and has the additional benefit of 
being simple, is all we really want for the study of finite nuclei. In the next section we show 
that in our model, a simple Skyrme-like functional satisfies these requirements. 



IV. PHENOMENOLOGICAL FUNCTIONAL AND SKYRME INTERACTIONS 



We want a Skyrme-like intrinsic-density functional with a minimum at the correct value 
to next-to- leading order in 1/N, that is, a functional that at its minimum reproduces the 
0(N 3 ) and 0(N 2 ) terms (the latter of which is zero) in Eq. (jUJ), and the second line of Eq. 
(|1(J|) when pi is calculated via Eq. (J15)) . We will reproduce the intrinsic density if the <fi that 
solves the Kohn-Sham equation is 

4>{x) = o (z) + ~M X ) ( 31 ) 



N 

y/{N-l)mg 
2 

y/(N-l)mg 



cosh 



.Nmgx 



-cosh 



2c N 
_ 3 Nmgx 



2c 



— cosh 



N 



, Nmgx , 
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where s = 1 + m 2 g 2 N 2 p gs /2 = 1 + vr 2 /6, and the cat's appear because of their presence in 
Eq. dHJ. 

To express such a functional through the formalism detailed above, we need to keep the 
center-of-mass kinetic energy in the Hamiltonian and use an effective nucleon mass 

m " = T^N > < 32 > 

where a is a constant to be determined later 3 . We denote by Eq[pi] the mean- field functional 
obtained by changing E Q [pj] in Eq. (fTHj) in that way (both changes affect the kinetic term), 
and write the functional through first order in 1/N as 

E[pj] = E*[ PI ] + EM , E l \p I ] = - fdx (||p 2 + ^-^P?) , (33) 

where (3 and 7 are also constants to be determined later. Eq. (|19j) . with the omission of all 
but the first three terms in Eq. (}2T)j) (which we argued could be neglected without changing 
the minimum of the functional), leads to the explicit Kohn- Sham-like equation 

-(N-l + f3)g\<P\ 2 <P-—L-\<P\ A <P = e<p. (34) 



2m* "' 2m* N 

Inserting <p = 4> + 4>i/N into this equation, writing e = e + ei/N, and expanding in inverse 
powers of (N — 1) rather than N (because of the c^r's in Eq. (|31j) h we find at leading order 
in (N — 1) _1 the self-consistent equation 

~<I>1 -(N- l)(7|0o| 2 0o = eo0o • (35) 
At 0(1/[N - 1]) = 0(1/N), assuming <p real, we have 

- 3(JV - l)^^0i - —</>» -(N- l)g^l - = e o 0i + (36) 

Now, the solution to Eq. §^ is the O of Eq. (jSH), with e = -mg 2 (N - l) 2 /8. The 
solution to Eq. (|3T)j). with the appropriate 4> and e plugged in, is the 0i of Eq. (|3*T|) . provided 
(3 = a — 2s, 7 = 18s, and e\ = ae . The requirement that the energy E[pi\ in Eq. (|3^|) have 
no 0(N 2 ) term at the minimum is satisfied as well if and only if 

t* = 2(l + y), /? = 2(1 - |) , 7 = 18s. (37) 

Thus, minimizing our functional in Eq. (J33j) with these values of the parameters yields the 
correct energy and intrinsic density to 0(1/N). 

The functional can be extended to reproduce pi and E at higher orders in (1/N) with 
the addition of higher powers of pi/N 2 , though finding the corresponding coefficients gets 
harder. These functionals, as discussed, do not converge to the "exact one" for all pi, but 
they do around the ground-state pi. The result suggests that real Skyrme functionals, which 



3 Though an effective mass can be included through the use of an intrinsic kinetic-energy density 13], we 
will simply alter the nucleon mass here. 
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may be viewed as depending on intrinsic semilocal kinetic, spin-orbit, and current densities 
in addition to the ordinary intrinsic density, behave the same way; we have no reason to 
think them correct away from the minimum. 

It's worth noting that the construction above is even more straightforward if we take 
Pi(x) = N\(p(x + x)\ 2 rather than define it as the expectation value of pi in the the Slater 
determinant constructed from 0, i.e., as in Eq. (|15|) . We do the latter only because the mean- 
field approximation and, more importantly, real Skyrme junctionals work better with that 
definition li]|. The success of the wave-function-based assignment pi = (p/) m f raises the 



hope that Skyrme functionals can be derived and improved through mean-field theory with 
some effective low-momentum interaction, with the addition of relatively small corrections. 
The density-matrix expansion works along those lines, and deserves more attention. It 
remains to be seen whether the formalism presented here — the inversion method modified 
for self-bound systems — can be happily married to such an approach. 

We have said nothing here about "symmetry restoration", a technique that, particularly 
in its rotational form, is viewed as a step beyond the local or semi-local Skyrme-mean- 
field equations. The projected-mean-field equations, while more complicated than their 
unprojected counterparts, are still equations for single-particle orbitals, however. Symmetry- 
restored mean-field theory is simply a more accurate starting point for the inversion method 
than plain mean-field theory. Presumably, additional corrections can be absorbed into those 
equations, as they were here. This statement, of course, says nothing about what the 
corrected functionals will look like. 

To summarize, we have shown that a version of the Hohenberg-Kohn theorem holds 
for the intrinsic density in self-bound system and that, to the extent that a convergent 
expansion around mean-field solutions exists for such systems, a Kohn- Sham-like procedure 
can be applied to obtain functionals of the intrinsic density that include all correlations. 
Skyrme functionals seem to be approximations to the exact intrinsic-density functionals but 
it is quite possible that, like our phenomenological functional, they work poorly away from 
the ground-state density and energy. On the other hand, we do not always need them to 
do more than that, and a concerted effort to improve even this limited kind of functional 
is worthwhile. The density-matrix expansion, which can clearly be applied more accurately 
than it was 35 years ago, may offer improvements and should be further explored. 

Fruitful discussions with A. Bhattacharyya, T. Duguet, R.J. Furnstahl, T. Schaefer, and 
A. Schwenk are gratefully acknowledged. This work was supported in part by the U.S. 
Department of Energy under Contract DE-FG02-97ER41019. 
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